Comparison of NPS WBM Products

Comparing NPS WBM Products

This RMarkdown explores the comparability of several data products associated with the NPS water balance model (WBM):

1. WBM outputs generated from the WBM R code (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw)
2. WBM outputs generated for park centroids (https://parkfutures.s3.us-west-2.amazonaws.com/park-centroid-wb/index.html)
3. WBM outputs generated across CONUS grid cells (http://www.yellowstone.solutions/thredds/catalog/daily_or_monthly/v2_historical/gridmet_v_1_5_historical/catalog.html)

We will be using Petersburg National Battlefield (PETE) as our location to compare each of these products’ outputs.

WBM R Code

For this comparison, we used the R scripts associated with the WBM found on AWS (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw). These folders contain code to calculate WBM variables for PETE. In the example provided, 10 random points within the park centroid’s MACA grid cell are selected, then used to run the WBM. At the end of the workflow, the 10 separate runs are averaged for each day.

First, let’s load in the functions associated with running the WBM in R:

source("aws_r_tutorial/wbm_functions/ET_functions.R")
source("aws_r_tutorial/wbm_functions/WB_functions.R")

Next, we can load in the metadata associated with the 10 randomly-selected points. This metadata is used as inputs in the WBM and is generated in a separate R script ahead of calculating the WBM.

wb_sites <- read.csv("aws_r_tutorial/PETE_site_parameters.csv") 

We also load in the historical GridMET variables associated with our centroid points in PETE. Because we are currently unable to locate the input dataset used in the tutorial hosted on AWS (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw), we are instead using the climate data hosted on AWS for the PETE’s centroid (https://parkfutures.s3.us-west-2.amazonaws.com/maca-tprh-data/index.html). Moreover, we are choosing to start executing the WBM in 1979; we are under the impression that the year 1979 is used as a “warm up” period for both the centroid and the gridded data sets.

Historical <- read.csv("aws_data_downloads/centroid_data/PETE_historical.csv")

Now, we are ready to run the WBM! The code below was pulled directly from the tutorial folders hosted on AWS (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw). In it, the Oudin method is used for calculating daily PET.

n <- nrow(wb_sites)
#Threshold temperature (deg C) for growing degree-days calculation
T.Base = 0 

#Method for PET calculation 
Method = "Oudin"  #Oudin is default method for daily PRISM and MACA data (containing only Tmax, Tmin, and Date). 

#Date format
DateFormat = "%m/%d/%Y"

############################################################ CREATE CLIMATE INPUTS #############################################################
#### Historical
# Convert pr.In to mm and F to C
Historical$ppt_mm <- (Historical$Precip..in.*25.4)
Historical$tmax_C <- 5/9*(Historical$Tmax..F. - 32)
Historical$tmin_C <- 5/9*(Historical$Tmin..F. - 32)
Historical$tmean_C <- 5/9*(Historical$Tavg..F. - 32)

######################################################### END CLIMATE INPUTS ####################################################################


######################################################### CALCULATE WB VARIABLES ################################################################
AllDailyWB<-list()
DailyWB<-Historical

for(i in 1:nrow(wb_sites)){
  ID = wb_sites$SiteID[i] # KW: was previously "Site"
  Lat = wb_sites$Lat[i]
  Lon = wb_sites$Lon[i]
  Elev = wb_sites$Elev[i] # KW: was previously "Elevation"
  Aspect = wb_sites$Aspect[i]
  Slope = wb_sites$Slope[i]
  SWC.Max = wb_sites$SWC.Max[i]
  Wind = wb_sites$Wind[i]
  Snowpack.Init = wb_sites$Snowpack.Init[i]
  Soil.Init = wb_sites$Soil.Init[i]
  Shade.Coeff = wb_sites$Shade.Coeff[i]
  
  #Calculate daily water balance variables 
  
  DailyWB$ID = ID
  DailyWB$doy <- yday(DailyWB$Date)
  DailyWB$daylength = get_daylength(DailyWB$Date, Lat)
  DailyWB$jtemp = as.numeric(get_jtemp(Lon, Lat))
  DailyWB$F = get_freeze(DailyWB$jtemp, DailyWB$tmean_C)
  DailyWB$RAIN = get_rain(DailyWB$ppt_mm, DailyWB$F)
  DailyWB$SNOW = get_snow(DailyWB$ppt_mm, DailyWB$F)
  DailyWB$MELT = get_melt(DailyWB$tmean_C, DailyWB$jtemp, hock=4, DailyWB$SNOW, Snowpack.Init)
  DailyWB$PACK = get_snowpack(DailyWB$jtemp, DailyWB$SNOW, DailyWB$MELT)
  DailyWB$W = DailyWB$MELT + DailyWB$RAIN
  if(Method == "Hamon"){
    DailyWB$PET = ET_Hamon_daily(DailyWB)
  } else {
    if(Method == "Penman-Monteith"){
      DailyWB$PET = ET_PenmanMonteith_daily(DailyWB)
    } else {
      if(Method == "Oudin"){
        DailyWB$PET = get_OudinPET(DailyWB$doy, Lat, DailyWB$PACK, DailyWB$tmean_C, Slope, Aspect, Shade.Coeff)
      } else {
        print("Error - PET method not found")
      }
    }
  }
  DailyWB$PET = modify_PET(DailyWB$PET, Slope, Aspect, Lat, Shade.Coeff)
  DailyWB$W_PET = DailyWB$W - DailyWB$PET
  DailyWB$SOIL = get_soil(DailyWB$W, Soil.Init, DailyWB$PET, DailyWB$W_PET, SWC.Max)
  DailyWB$DSOIL = diff(c(Soil.Init, DailyWB$SOIL))
  DailyWB$AET = get_AET(DailyWB$W, DailyWB$PET, DailyWB$SOIL, Soil.Init)
  DailyWB$W_ET_DSOIL = DailyWB$W - DailyWB$AET - DailyWB$DSOIL
  DailyWB$D = DailyWB$PET - DailyWB$AET
  DailyWB$GDD = get_GDD(DailyWB$tmean_C, T.Base)
  AllDailyWB[[i]] = DailyWB
}

WBData<-do.call(rbind,AllDailyWB)

Lastly, we aggregate the 10 WBM outputs by date and grab the means for a handful of the WBM variables we want to compare. We also convert the output (mm) to inches.

WBData_coded <- do.call(rbind, AllDailyWB) %>%
  group_by(Date, GCM) %>%
    # convert to inches
  dplyr::summarize(mean_runoff = mean(W_ET_DSOIL, na.rm = TRUE)/25.4,
                   mean_pet = mean(PET, na.rm = TRUE)/25.4,
                   mean_rain = mean(RAIN, na.rm = TRUE)/25.4,
                   mean_snow = mean(SNOW, na.rm = TRUE)/25.4,
                   mean_melt = mean(MELT, na.rm = TRUE)/25.4,
                   mean_snowpack = mean(PACK, na.rm = TRUE)/25.4,
                   mean_aet = mean(AET, na.rm = TRUE)/25.4) %>%
  # Remove "warm up" year 
  filter(year(as_date(Date)) > 1979)

Next, we can pull in and organize the data associated with the other data products on AWS: the gridded product, and the centroid product.

Park Centroid Data

Here is the centroid data downloaded directly from AWS for PETE (https://parkfutures.s3.us-west-2.amazonaws.com/park-centroid-wb/index.html). We also remove the first year of data (1979), since we are under the impression that this year is used as the “warm up” period, and because our gridded data set only has data starting in 1980. Data are reported in inches in this data set. We are also under the impression that the Oudin method was selected for calculating PET.

WBData_online_centroid <- read_csv("aws_data_downloads/centroid_data/PETE_water_balance_historical.csv") %>%
dplyr::filter(year(as_date(Date)) > 1979)

Gridded CONUS Data

Here, we pull in the gridded data, which was downloaded directly from AWS, then clipped to PETE’s park boundary. Gridded data starts in 1980 (which could be explained by the fact that 1979 is used as a “warm up”). We are also under the impression that the Oudin method was selected for calculating PET. The python script provided to us that represents the workflow used to develop the gridded dataset can be found at “aws_python_script/tercek_script.py”.

WBData_online_grid <- rast("aws_data_downloads/gridded_data/PETE_runoff_historic_gridmet_1980_2023.tif")

# next we clip the grid to only grids that intersect our wb_sites points in PETE:
NPS_points <- st_transform(sf::st_as_sf(wb_sites, coords = c("Lon", "Lat"), crs = 4326), st_crs(WBData_online_grid))
grid_runoff <- terra::extract(WBData_online_grid, NPS_points) %>%
  pivot_longer(cols = -ID) %>%
  dplyr::group_by(name) %>%
    # convert to inches
  dplyr::summarize(mean_runoff = mean(value, na.rm = TRUE)/25.4/10) %>%
  filter(name != "runoff_NA") %>%
  mutate(date = as_date(sub(".*_", "", name))) %>%
  filter(year(date) < 2023)

WBData_online_grid <- rast("aws_data_downloads/gridded_data/PETE_pet_historic_gridmet_1980_2023.tif")  
  
grid_pet <- terra::extract(WBData_online_grid, NPS_points) %>%
  pivot_longer(cols = -ID) %>%
  dplyr::group_by(name) %>%
  # convert to inches
  dplyr::summarize(mean_pet = mean(value, na.rm = TRUE)/25.4/10) %>%
  filter(name != "PET_NA") %>%
  mutate(date = as_date(sub(".*_", "", name))) %>%
  filter(year(date) < 2023)

WBData_online_grid <- rast("aws_data_downloads/gridded_data/PETE_rain_historic_gridmet_1980_2023.tif")

grid_rain <- terra::extract(WBData_online_grid, NPS_points) %>%
  pivot_longer(cols = -ID) %>%
  dplyr::group_by(name) %>%
    # convert to inches
  dplyr::summarize(mean_rain = mean(value, na.rm = TRUE)/25.4/10) %>%
  filter(name != "rain_NA") %>%
  mutate(date = as_date(sub(".*_", "", name))) %>%
  filter(year(date) < 2023)

Comparing Outputs

Runoff

Here I am plotting each of the different WBM runoff products against each other for comparison.

year_colors <- c("#002EA3", "#C3F73A", "#E70870", "#256BF5", "#745CFB", "#FFCA3A", "#578010", "#56104E")
interpolated_colors <- colorRampPalette(year_colors)(length(unique(year(WBData_coded$Date))))

  one <- ggplot() +
    ggtitle("Runoff in") +
    theme_bw() +
    xlab("Generated Data from Code") +
    ylab("Centroid Data on AWS") +
    geom_point(aes(x = WBData_coded$mean_runoff, 
                   y = WBData_online_centroid$`runoff in`, 
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

 two <-  ggplot() +
    theme_bw() +
    xlab("Generated Data from Code") +
    ylab("Gridded Data on AWS") +
    geom_point(aes(x = WBData_coded$mean_runoff,
                   y = grid_runoff$mean_runoff,
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

 three <-  ggplot() +
    theme_bw() +
    xlab("Centroid Data on AWS") +
    ylab("Gridded Data on AWS") +
    geom_point(aes(x = WBData_online_centroid$`runoff in`,
                   y = grid_runoff$mean_runoff,
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

ggpubr::ggarrange(one + theme(legend.position = "none"),
                  two + theme(legend.position = "none"),
                  three + theme(legend.position = "none"), 
                  ncol = 1, nrow = 3,
  common.legend = TRUE, legend = "right")

All three runoff products plotted together through time:

ggplotly(
  ggplot() +
    geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `runoff in`), color = "#E70870", alpha = 1) +
    geom_point(data = grid_runoff, aes(x = date, y = mean_runoff), color = "#256BF5", alpha = 0.6) +
    geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_runoff), color = "black", cex = 0.2) +
    theme_bw() +
    labs(x = "Date", y = "Runoff in") 
    # geom_text(data = WBData_online_centroid[1, ], aes(x = as.Date(Date), y = `runoff in`), label = "Centroid Data on AWS", color = "#E70870", vjust = -35, hjust = -1) +
    # geom_text(data = grid_runoff[1, ], aes(x = date, y = mean_runoff), label = "Gridded Data on AWS", color = "#256BF5", vjust = -33, hjust = -1) +
    # geom_text(data = WBData_coded[1, ], aes(x = as.Date(Date), y = mean_runoff), label = "Generated Data from Code", color = "black", vjust = -31, hjust = -0.8)
    )

Figure Legend: Centroid Data on AWS, Generated Data from Code, Gridded Data o AWS.

PET

Here I am plotting each of the different WBM PET products against each other for comparison.

 one <- ggplot() +
    ggtitle("PET in") +
    theme_bw() +
    xlab("Generated Data from Code") +
    ylab("Centroid Data on AWS") +
    geom_point(aes(x = WBData_coded$mean_pet, 
                   y = WBData_online_centroid$`PET in`, 
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

 two <-  ggplot() +
    theme_bw() +
    xlab("Generated Data from Code") +
    ylab("Gridded Data on AWS") +
    geom_point(aes(x = WBData_coded$mean_pet,
                   y = grid_pet$mean_pet,
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

 three <-  ggplot() +
    theme_bw() +
    xlab("Centroid Data on AWS") +
    ylab("Gridded Data on AWS") +
    geom_point(aes(x = WBData_online_centroid$`PET in`,
                   y = grid_pet$mean_pet,
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

ggpubr::ggarrange(one + theme(legend.position = "none"),
                  two + theme(legend.position = "none"),
                  three + theme(legend.position = "none"),  
                  ncol = 1, nrow = 3,
  common.legend = TRUE, legend = "right")

All three PET products plotted together through time:

ggplotly(
  ggplot() +
    geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `PET in`), color = "#E70870", alpha = 1) +
    geom_point(data = grid_pet, aes(x = date, y = mean_pet), color = "#256BF5", alpha = 0.6) +
    geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_pet), color = "black", cex = 0.2) +
    theme_bw() +
    labs(x = "Date", y = "PET in")
)

Figure Legend: Centroid Data on AWS, Generated Data from Code, Gridded Data on AWS.

Rain

Here I am plotting each of the different WBM rain products against each other for comparison.

 one <- ggplot() +
    ggtitle("Rain in") +
    theme_bw() +
    xlab("Generated Data from Code") +
    ylab("Centroid Data on AWS") +
    geom_point(aes(x = WBData_coded$mean_rain, 
                   y = WBData_online_centroid$`rain in`, 
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

 two <-  ggplot() +
    theme_bw() +
    xlab("Generated Data from Code") +
    ylab("Gridded Data on AWS") +
    geom_point(aes(x = WBData_coded$mean_rain,
                   y = grid_rain$mean_rain,
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

 three <-  ggplot() +
    theme_bw() +
    xlab("Centroid Data on AWS") +
    ylab("Gridded Data on AWS") +
    geom_point(aes(x = WBData_online_centroid$`rain in`,
                   y = grid_rain$mean_rain,
                   color = factor(year(WBData_coded$Date))),
               size = 1) +
    scale_color_manual(values = interpolated_colors) +
    labs(color = "Year")

ggpubr::ggarrange(one + theme(legend.position = "none"), 
                  two + theme(legend.position = "none"), 
                  three + theme(legend.position = "none"),  
                  ncol = 1, nrow = 3,
  common.legend = TRUE, legend = "right")

All three rain products plotted together through time:

  ggplotly(
    ggplot() +
    geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `rain in`), color = "#E70870", alpha = 1) +
    geom_point(data = grid_rain, aes(x = date, y = mean_rain), color = "#256BF5", alpha = 0.6) +
    geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_rain), color = "black", cex = 0.2) +
    theme_bw() +
    labs(x = "Date", y = "Rain in")
  )

Figure Legend: Centroid Data on AWS, Generated Data from Code, Gridded Data o AWS.

Discrepancies in outputs

Based on discussions with the NPS group, we understand that the published centroid data and the gridded data products were developed with different intentions, and subsequently, are different by design. For example:

• The gridded dataset’s workfow has the soil water storage start at full capacity, whereas the WBM R code has the soil water storage start at 0.

• The WBM R code’s calculation for Oudin PET is slightly different: it incorporates a shade coefficient to reduce PET (though it is set to NULL and therefore not in use), whereas the gridded workflow does not include this option. The WBM R code also sets PET to 0 if there is at least 2 mm of snowpack, whereas the gridded product sets PET to 0 if there is any snow (snowpack > 0).

However, we are under the impression that the park centroid data was generated using the same workflow outlined in the WBM R code. Though they are incredibly similar, they are not identical. Here, we lay out a few potential explanations for the differences in their outputs.

1. In the WBM code, we found a potential redundancy in calculating PET. It appears that get_OudinPET() incorporates aspects of modify_PET(). Both are used in the published tutorial on AWS in generating WBM variables. (See lines 107-120 above.) Below are the two functions:

get_OudinPET = function(doy, lat, snowpack, tmean, slope, aspect, shade.coeff=NULL){
  d.r = 1 + 0.033*cos((2*pi/365)*doy)
  declin = 0.409*sin((((2*pi)/365)*doy)-1.39)
  lat.rad = (pi/180)*lat
  sunset.ang = acos(-tan(lat.rad)*tan(declin))
  R.a = ((24*60)/pi)*0.082*d.r*((sunset.ang*sin(lat.rad)*sin(declin)) + (cos(lat.rad)*cos(declin)*sin(sunset.ang)))
  Oudin = ifelse(snowpack>2,0,ifelse(tmean>-5,(R.a*(tmean+5)*0.408)/100,0))
  Folded_aspect = abs(180-abs((aspect)-225))
  Heatload = (0.339+0.808*cos(lat*(pi/180))*cos(slope*(pi/180)))-(0.196*sin(lat.rad)*sin(slope*(pi/180)))-(0.482*cos(Folded_aspect*(pi/180))*sin(slope*(pi/180)))
  sc = ifelse(!is.null(shade.coeff), shade.coeff, 1)
  OudinPET = Oudin * Heatload * sc
  return(OudinPET)
}

modify_PET = function(pet, slope, aspect, lat, freeze, shade.coeff=NULL){
  f.aspect = abs(180 - abs(aspect - 225))
  lat.rad = ifelse(lat > 66.7, (66.7/180)*pi, (lat/180)*pi)
  slope.rad = (slope/180)*pi
  aspect.rad = (f.aspect/180)*pi
  heat.load = 0.339+0.808*cos(lat.rad)*cos(slope.rad) - 0.196*sin(lat.rad)*sin(slope.rad) - 0.482*cos(aspect.rad)*sin(slope.rad)
  sc = ifelse(!is.null(shade.coeff), shade.coeff, 1)
  freeze = ifelse(freeze == 0,0,1)
  PET.Lutz = pet*heat.load*sc*freeze
  return(PET.Lutz)
}

Of note, the exclusion of the modify_PET() step does not lead to a major change in values. Below is a plot of the difference between a runoff output that uses modify_PET() and the runoff output that does not use modify_PET():

n <- nrow(wb_sites)
#Threshold temperature (deg C) for growing degree-days calculation
T.Base = 0 

#Method for PET calculation 
Method = "Oudin"  #Oudin is default method for daily PRISM and MACA data (containing only Tmax, Tmin, and Date). 

#Date format
DateFormat = "%m/%d/%Y"

############################################################ CREATE CLIMATE INPUTS #############################################################
#### Historical
# Convert pr.In to mm and F to C
Historical$ppt_mm <- (Historical$Precip..in.*25.4)
Historical$tmax_C <- 5/9*(Historical$Tmax..F. - 32)
Historical$tmin_C <- 5/9*(Historical$Tmin..F. - 32)
Historical$tmean_C <- 5/9*(Historical$Tavg..F. - 32)

######################################################### END CLIMATE INPUTS ####################################################################


######################################################### CALCULATE WB VARIABLES ################################################################
AllDailyWB<-list()
DailyWB<-Historical

for(i in 1:nrow(wb_sites)){
  ID = wb_sites$SiteID[i] # KW: was previously "Site"
  Lat = wb_sites$Lat[i]
  Lon = wb_sites$Lon[i]
  Elev = wb_sites$Elev[i] # KW: was previously "Elevation"
  Aspect = wb_sites$Aspect[i]
  Slope = wb_sites$Slope[i]
  SWC.Max = wb_sites$SWC.Max[i]
  Wind = wb_sites$Wind[i]
  Snowpack.Init = wb_sites$Snowpack.Init[i]
  Soil.Init = wb_sites$Soil.Init[i]
  Shade.Coeff = wb_sites$Shade.Coeff[i]
  
  #Calculate daily water balance variables 
  
  DailyWB$ID = ID
  DailyWB$doy <- yday(DailyWB$Date)
  DailyWB$daylength = get_daylength(DailyWB$Date, Lat)
  DailyWB$jtemp = as.numeric(get_jtemp(Lon, Lat))
  DailyWB$F = get_freeze(DailyWB$jtemp, DailyWB$tmean_C)
  DailyWB$RAIN = get_rain(DailyWB$ppt_mm, DailyWB$F)
  DailyWB$SNOW = get_snow(DailyWB$ppt_mm, DailyWB$F)
  DailyWB$MELT = get_melt(DailyWB$tmean_C, DailyWB$jtemp, hock=4, DailyWB$SNOW, Snowpack.Init)
  DailyWB$PACK = get_snowpack(DailyWB$jtemp, DailyWB$SNOW, DailyWB$MELT)
  DailyWB$W = DailyWB$MELT + DailyWB$RAIN
  if(Method == "Hamon"){
    DailyWB$PET = ET_Hamon_daily(DailyWB)
  } else {
    if(Method == "Penman-Monteith"){
      DailyWB$PET = ET_PenmanMonteith_daily(DailyWB)
    } else {
      if(Method == "Oudin"){
        DailyWB$PET = get_OudinPET(DailyWB$doy, Lat, DailyWB$PACK, DailyWB$tmean_C, Slope, Aspect, Shade.Coeff)
      } else {
        print("Error - PET method not found")
      }
    }
  }
  #DailyWB$PET = modify_PET(DailyWB$PET, Slope, Aspect, Lat, Shade.Coeff)
  DailyWB$W_PET = DailyWB$W - DailyWB$PET
  DailyWB$SOIL = get_soil(DailyWB$W, Soil.Init, DailyWB$PET, DailyWB$W_PET, SWC.Max)
  DailyWB$DSOIL = diff(c(Soil.Init, DailyWB$SOIL))
  DailyWB$AET = get_AET(DailyWB$W, DailyWB$PET, DailyWB$SOIL, Soil.Init)
  DailyWB$W_ET_DSOIL = DailyWB$W - DailyWB$AET - DailyWB$DSOIL
  DailyWB$D = DailyWB$PET - DailyWB$AET
  DailyWB$GDD = get_GDD(DailyWB$tmean_C, T.Base)
  AllDailyWB[[i]] = DailyWB
}

WBData_coded_no_modify_PET <- do.call(rbind, AllDailyWB) %>%
  group_by(Date, GCM) %>%
    # convert to inches
  dplyr::summarize(mean_runoff = mean(W_ET_DSOIL, na.rm = TRUE)/25.4,
                   mean_pet = mean(PET, na.rm = TRUE)/25.4,
                   mean_rain = mean(RAIN, na.rm = TRUE)/25.4,
                   mean_snow = mean(SNOW, na.rm = TRUE)/25.4,
                   mean_melt = mean(MELT, na.rm = TRUE)/25.4,
                   mean_snowpack = mean(PACK, na.rm = TRUE)/25.4,
                   mean_aet = mean(AET, na.rm = TRUE)/25.4) %>%
  # Remove "warm up" year 
  filter(year(as_date(Date)) > 1979)

ggplot() +
  #geom_point(data = WBData_coded, aes(x = as_date(Date), y = mean_runoff), color = "black", cex = 1.5) +
  geom_col(aes(x = as_date(WBData_coded$Date), y = WBData_coded$mean_runoff - WBData_coded_no_modify_PET$mean_runoff), color = "#578010") +
  theme_bw() +
  ggtitle("Difference in Runoff (inches)") +
  xlab("Date") +
  ylab("With - Without `modify_PET()`")

2. Though we assumed that the coordinates provided in the folder hosted on AWS (“aws_r_tutorial/PETE_site_parameters.csv”) were used in developing PETE’s published centroid data, it is possible that this is not the case.

3. We were unable to find the exact historical climate dataset used in the AWS tutorial. Instead, we used PETE’s centroid’s historical GridMET data posted on AWS (https://parkfutures.s3.us-west-2.amazonaws.com/maca-tprh-data/index.html). If a different data set was used to generate the centroid data, this could also explain the discrepancies.

4. The code we are running starts the WBM in the year 1979. Though we believe that 1979 is also the year used to generate the centroid data, it is possible another year was used as the starting point.